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Abstract. During the recovery phase of a magnetic storm, fluxes of relativistic 
(> 1 MeV) electrons in the inner magnetosphere (3 < L < 6) increase to beyond 
pre-storm levels, reaching a peak about 4 days after the initiation of the storm. In order 
to account for the generation of these "killer electrons" , a model is presented primarily 
based on stochastic acceleration of electrons by enhanced whistler-mode chorus. In 
terms of a quasi-linear formulation, a kinetic (Fokker-Planck) equation for the electron 
energy distribution is derived, comprising an energy diffusion coefficient based on 
gyroresonant electron-whistler-mode wave interaction and parallel wave propagation; a 
source term representing substorm-produced (lower energy) seed electrons; and a loss 
term representing electron precipitation due to pitch-angle scattering by whistler-mode 
waves and EMIC waves. Steady-state solutions for the electron energy distribution are 
constructed, and fitted to an empirically-derived relativistic Maxwellian distribution 
for the high energy "hard" electron population at geosynchronous orbit. If the average 
whistler amplitude is sufficiently large, for instance 75 pT - 400 pT, dependent on the 
values of the other model parameters, and assuming a background plasma density of 
A^o = 10 cm~^ outside the plasmasphere, then a good fit to the empirical distribution 
is obtained, and corresponds to a timescale for the formation of the high-energy steady 
state distribution of 3 - 5 days. For a lower representative value of the background 
plasma density, A'^q = 1 cm~^, smaller whistler amplitudes, in the range 13 - 72 pT, 
can produce the high-energy distribution in the required time frame of several days. It 
is concluded from the model calculations that the process of stochastic acceleration by 
gyroresonant electron- whistler-mode wave interaction, in conjunction with pitch- angle 
scattering by EMIC waves, constitutes a viable mechanism for generating "killer 
electrons" during geomagnetic storms. The mechanism is expected to be particularly 
effective for the class of small and moderate storms possessing a long-lasting recovery 
phase during which many substorms occur. 



1. INTRODUCTION 

It is well known that variations in the fluxes of relativistic electrons, of kinetic 
energies > 1 MeV, in the inner magnetosphere (3 < L < 6) are related to disturbed 
magnetospheric conditions commonly called "magnetic storms". Typically, for many 
storms the electron fluxes diminish rapidly during the main phase of the storm . The 
main phase depletion of relativistic electrons occurs in association with large negative 
values of the interplanetary magnetic field B^ and large sudden increases in the solar 
wind density and pressure \Faulikas and Blake, 1979| ; \Blake et ai, 1997| . Subsequently, 



during the recovery phase of the storm, fluxes increase to beyond pre-storm levels and 
peak about 4 days after the initiation of the storm \\Paulikas and Blake, 1979|; \Bake7 



\et ai, 1986, 1994a, 199^ ; \^agai, 1988t \Li et ai, 1997at \Keeves et ai, 1998|| . These 



enhancements in fluxes of relativistic electrons, which are colloquially referred to as 
"killer electrons" , have become the subject of considerable attention by magnetospheric 
physicists. Not only do the enhancements constitute an intrinsically interesting physics 
problem in the near-Earth space, but they constitute a potentially serious hazard to 
satellites, space stations, and, conceivably, humans in space. In fact, satellite disfunctions 
( "anomalies" ) have been linked to the effects of relativistic electron increases Waker ei 



\al., 1994b, 1997|| , and the state of the radiation belt environment has become a major 
concern in space weather forecasting [e.g.. Baker, 1998; Reeves, 1998a]. The region near 
geosynchronous (or geostationary) orbit, L ~ 6.6, in the geographic equatorial plane, 
is of particular interest because it is the operating zone of many orbiting satellites. 
Reeves [1998b] has recently examined the relationship between relativistic electron 
enhancements at geosynchronous orbit and magnetic storms as measured by the Dgt 
index. In particular, the 30 most intense relativistic electron events during 1992-1995 
were examined, and it was found that every relativistic electron event was associated 
with a magnetic storm as indicated by the Dgt index, though a small fraction (about 
10%) of magnetic storms did occur with no increase in relativistic electron fluxes. Thus, 



one conclusion from Reeves' [1998b] analysis is that intense solar wind conditions are 
necessary to generate strong relativistic electron enhancements. Nevertheless, despite the 
accumulated magnetic storm data from satellites over many years, including coordinated 
observations from the International Solar- Terrestrial Physics (ISTP) constellation of 
spacecraft and other multi-satellite missions [e.g.. Baker et ai, 1997; Reeves et al., 
1998], there is, as yet, no accepted explanation for the generation of the relativistic 
electrons. Specifically, it is not known exactly how, where, or when the electrons are 
accelerated. Various energization mechanisms have been proposed, and most of these 
are reviewed by Li et al. [1997a]. It appears easier to explain the main phase depletion 
of energetic electrons than their subsequent recovery and enhancement. The drop in 
relativistic electron fluxes near geosynchronous orbit is partly due to adiabatic responses 
(conserving all three adiabatic invariants) to magnetic field decreases, as refiected in 
the reduction in Dst index [e.g., Kim and Chan., 1997]. Nevertheless, Li et al. [1997a] 
show that other physical mechanisms, including precipitation, must also contribute to 



the depletion. It has been suggested that radial diffusion \Schulz and Lanzerotti, 1974 
invoked to explain the existence of the outer electron radiation belt itself, could also 
generate electrons of MeV energies in the inner magnetosphere. This mechanism, which 
involves inwardly transporting energetic electrons from a presumed source in the outer 
magnetosphere (in the tail), can produce such energies during relatively quiet periods 
[e.g., Selesnick and Blake, 1997a], although the process is too slow during active times 
Li et al., 1997a| ; \Blake et al., 199^] . Certain global recirculation processes, involving 



radial diffusion, have also been proposed to generate relativistic electrons [e.g.. Baker et 
al, 1986, 1989; Fujimoto and Nishida, 1990], though these have proved inadequate, as 
the transport rates are too slow. Sheldon et al. [1998] have recently identified the outer 
polar cusp region as a potential acceleration region of the magnetosphere and possible 
source of energetic electrons for the outer radiation belt, though further calculations 
are needed to evaluate the significance of the study. In another mechanism still to 



be fully evaluated, Rostoker et al. [1998] and Liu et al. [1999] make the case that 
large-amplitude ULF pulsations have the potential to supply the energy necessary to 
create the enhanced relativistic electron fluxes. In a 3D global MHD simulation of the 
rapid rise of relativistic electron fluxes during the January 1997 magnetic cloud event, 
Hudson et al. [1999a, b] also found that ULF oscillations may play a role in energizing 
relativistic electrons, via a mechanism involving drift-resonant acceleration and radial 
transport. 

It is becoming increasingly apparent that electrons are accelerated to relativistic (> 
1 MeV) energies in situ in the inner magnetosphere [e.g., Blake et al., 1998]. Significant 
evidence in support of this conclusion is the observation by Selesnick and Blake [1997b] 
that the phase space density of electrons of greater than MeV energies peaks near L=4 
to L=5 during storms. As a result of substorm activity, electrons with energies up to 
~ 300 keV are injected near geosynchronous orbit \\Cayton et al., 1989| ; \Baker et al. 



1989, 1998| . These electrons appear to form the source population for the relativistic 
electrons of greater than MeV energies that are subsequently observed. Summers et al. 
[1998, 1999] have shown that whistler-mode waves could provide an effective mechanism 
for accelerating electrons from energies near 100 keV to above 1 MeV in the region 
outside the plasmapause during the storm recovery phase. In a survey of potential 
wave modes for electron scattering and stochastic acceleration to relativistic energies 
during magnetic storms. Home and Thome [1998] concluded, in particular, that in 
low density regions of the magnetosphere where the electron gyrofrequency exceeds 
the electron plasma frequency, there are four potential wave modes that can resonate 
with electrons in the energy range 100 keV to a few MeV: the whistler, LO, RX, and Z 
modes. The concept of stochastic acceleration of electrons by whistler-mode waves in 
the magnetosphere has also been discussed by Temerin et al. [1994], Li et al. [1997a], 
Temerin [1998], and Roth et al. [1999]. It is the purpose of present paper to quantify the 
model presented qualitatively in Section 8 of Summers et al. [1998] for the stochastic 
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acceleration of relativistic electrons during geomagnetic storms. Essential ingredients 
in the model are the spatial regions within the inner magnetosphere 3 < L < 9 where 
enhanced whistler-mode chorus \ Tsurutam and Smith, 1974, 1977| ; Koons and Roeder, 



19901 ; \Parrot and Gaye, 1994\\ and enhanced electromagnetic ion cyclotron (EMIC or 



L-mode) waves occur \\CornwaU et al., 19701 ; \Perraut et ai, 1976| ; \Jordanova et a/.] 



|1997| ; \Kozyra et al, 1997| ; see Figure 1. The aforementioned substorm-produced seed 



population of electrons in the energy range 100 keV to 300 keV have approximately 
circular drift paths within the region 3 < L < 9, and consequently will traverse the 
regions of enhanced whistler-mode chorus and enhanced EMIC waves. Specifically, 
in this paper we shall model the acceleration of electrons during the storm recovery 
phase by means of second-order Fermi (or stochastic gyroresonant) acceleration by weak 
whistler-mode turbulence. In a standard quasi-linear formulation, we construct a kinetic 
equation for the evolution of the electron energy distribution function. The equation 
contains an energy diffusion coefficient due to resonant whistler-mode wave/electron 
interaction, and an electron loss term due to pitch-angle scattering by both the 
whistler-mode and EMIC waves. We should point out here that the data from the Solar, 
Anomalous, and Magnetospheric Particle Explorer (SAMPEX) satellite \Nakam,ura et 



L, 1995| ; \Li et ai, 1997a| ; \Nakaniura, 1998| show that bursty electron precipitation 



occurs as the electron flux increases during the storm recovery phase. Such precipitation, 
together with the existence of an abundant supply of storm-produced lower-energy seed 
electrons \\Li et al, 1997a| \Baker et ai, 199"8| are supportive of the model constructed 



in this paper. The model is presented in detail in Section 2, and numerical solutions are 
presented in Section 3. The solutions are compared with data on the electron energy 
distribution of high energy (300 - 2000 keV) electrons at geosynchronous orbit. We find 
that stochastic gyroresonant acceleration by whistler-mode waves can indeed accelerate 
substorm-produced seed electrons in the inner magnetosphere to generate high-energy 
electron spectra of the type observed, following continuous injection of seed electrons 
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over a timescale of several days. Specific predictions of the model depend, of course, on 
the values taken for the model parameters. In Section 4 we briefiy assess our findings 
and state our conclusions. 



2. MODEL 

The region to which the model constructed in this paper applies is that part of 
the inner magnetosphere during storm-time that is illustrated in the idealized Figure 
|I|, for 3 < L < 9. This region contains an extensive subregion of whistler-mode 
chorus, a smaller but intense region of EMIC waves, and also contains the important 
geosynchronous-orbit region near L=6.6. Figure 1 is a simplified version of Figure 7 
of Summers et al. [1998], which itself was constructed on the basis of observations 
and relevant theory \Cornwall et al., 19701 ; \Perraut et ai, 1976| ; [Koons and Roeder, 



1990 ; Parrot and Gaye, 1994 ; Tsurutani and Smith, 1974, 1977 ; Kozyra et al.,\ 



19971 ; \Jordanova et al, 1997|] . We assume that geomagnetic storm activity produces a 



seed population of electrons of energy ?^ 100 keV as a source for the region specified 
in Figure 1. We are not concerned in this paper with the precise means (transport, 
original source location, etc.) by which the source is supplied. According to standard 
particle drift theory [e.g.. Wolf, 1995], the drift motion of "hot" (e.g., 100 keV) electrons 
close to the Earth is dominated by gradient drift, with the result that electrons 
execute approximately circular drift trajectories eastward about the Earth. Thus, the 
storm-supplied source electrons constitute a quasi-trapped population traversing the 
whistler-mode and EMIC wave subregions illustrated in Figure |I]. While executing the 
eastward drift, the electrons gyrate about the field lines and 'bounce' between mirror 
points, during which time they also undergo both energy and pitch-angle diffusion, as 
a result of their interaction with the whistler-mode chorus and EMIC waves. We shall 
assume that the pitch-angle scattering rate is much greater than the energy diffusion 
rate (see Appendix A), so that the electron distribution will be nearly isotropic. Further, 



Figure |I] 



we shall take account of pitch-angle diffusion of electrons into the loss-cone and their 
subsequent precipitation into the atmosphere, by characterizing their loss from the 
system by an escape time Tesc- The kinetic or Fokker-Planck equation describing the 
evolution in time t of the electron energy distribution function f{E, t) can be written as 

f = |l WE)f) - A [(^ME)- I ^. I) /] - £ + S(E, t), (1) 

where E = Ek/{meC^) = 7 — 1 is the particle kinetic energy in units of the rest mass 
energy, 7 = (1 — v'^ /c^)^^/"^ is the Lorentz factor, v is the particle speed, rrie is the 
electron rest mass, and c is the speed of light; f{E,t)dE is the number of particles per 
unit volume in the interval dE; D{E) is the energy diffusion coefficient due to resonant 
interaction of the electrons with whistler-mode turbulence; A[E) is the systematic 
acceleration rate due to the whistler-mode turbulence; | E^ \ is the energy loss rate 
due to processes not directly related to stochastic acceleration, namely, here assumed 
to be Coulomb collisions and synchrotron radiation; Tesc is the mean escape time of 
particles out of the system due to pitch-angle scattering by both whistler-mode and 
EMIC waves; and the source term S{E, t) represents the rate of particle injection into 
the inner magnetospheric region specified in Figure |l], as a result of storm activity. 
Equation (|l]) is not the standard form of Fokker-Planck equation employed in space 
physics, and so we give a brief account of its derivation in the Appendix. Detailed data 
on the whistler- mode chorus during storm-time are unfortunately not available. In fact, 
insufficient information is known about the energy spectrum of the turbulence in many 
space physics situations. While whistler-mode chorus emissions are normally considered 
to be discrete during geomagnetically quiet times \\Anderson and Kurth, 19"89|] , it is 
here assumed that during geomagnetic storms the concomitant, enhanced whistler-mode 
turbulence can be considered quasi-continuous. Specifically, a simplifying assumption is 
made that the whistler-mode turbulence is isotropic, homogeneous, stationary, and has 
a power-law spectral energy density distribution in wavenumber k, with spectral index 



g; specifically, the spectral energy density is assumed to take the form, 

W{k) = ^—^ (tn^Yw^ot, W,ot = r W{k)dk, (2) 

for wavenumbers greater than kmm-, to be specified below. In accordance with the 
quasi-linear diffusion formulation adopted in this paper, the whistler-mode turbulence 
is weak, i.e., comprises small-amplitude magnetic and electric wave fields. Momentum 
diffusion coefficients corresponding to whistler-mode waves have been obtained by 
various authors. We calculate the Fokker-Planck coefficients D{E) and A{E) from (A8) 
using the whistler-mode diffusion coefficients Dp derived by Hamilton and Petrosian 
[1992] (for 2 < g < 4), and Schlickeiser [1997] (for 1 < g < 2) for parallel wave 
propagation. The results are 

D{E) = Do [E{E + 2)f-'^/^ {E + l)-\ (3) 

A{E) = Doq[E{E + 2)f'^-'^/\ (4) 



where 



for 2 < g < 4; 



where 



7r(g 1) ( ckmin\ 2T3r, (k\ 



D{E)=Vo[E{E + 2)f\E + ir\ (6) 

A{E)=2Vo[E{E + 2)]-'/\ (7) 



8 V it / \m„ 



for 1 < g < 2. 

In (3) - (8), the two dimensionless parameters R and a are introduced; R is the ratio of 

turbulent energy Wtot to magnetic field energy, 

R = 87rWtot/B', = {AB/Bof, (9) 
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and 

a = nl/uj% = {mjm,)(3% (10) 

where fig = eBQ/{mf.c) is the electron gyrofrequency, with Bq the ambient magnetic 
field strength and e the electron charge; AB is the average whistler-mode wave 
amplitude; tUpe = (AirNQe^ /me)^^'^ is the electron plasma frequency, with A^o the 
particle number density; Jw is a weakly varying function of E; and j3a = Va/c where 
Va = Bq / (An NqUIpY^'^ is the Alfven speed, with nip the proton rest mass. The parameter 
a defined in (10) is identical to the parameter a used by Summers et al. [1998]. We 
note, in particular, as should indeed be the case, that since the kinetic energy variable 
E is dimensionless, the dimension of the diffusion coefficient D equals the dimension of 
the parameter Dq (or T?q) equals [time]^^. For definiteness, in (5) and (8) we set 

kmin = ^p/{cPa), (11) 

[e.g., Hamilton and Petrosian, 1992], where Qp = eBo/lmpc) is the proton gyrofrequency. 
The energy loss term | E^ |, in which we include losses due to Coulomb collisions 
and synchrotron radiation, can be expressed in the form 

\El\=6x IQ-^^NoiE + l)[E{E + 2)]"^/^ ^ ^33 x 10-^B^E{E + 2). (12) 

The first term on the right-hand side of the equation (12) is the energy loss rate due 
to Coulomb collisions, given by Melrose [1980], and the second term the energy loss 
rate due to synchrotron radiation, given by Blumenthal and Gould [1970]. In (12), the 
particle number density A'^o is in cm^'^, the magnetic field strength Bq is in gauss, and 
I El I is in sec""*". 

We note that there are potentially four infiuential parameters in the model: the 
parameter a defined by (10), the spectral index g, the turbulent wave power parameter 
R, and the mean particle escape-time Tesc- The value of the parameter a depends on 
the values taken for the particle density A'"o and the ambient magnetic field strength 
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Bq. We shall discuss Tesc, which we regard as an adjustable parameter, and the particle 
source function S below. The diffusion parameters Dq and "Do occurring in expressions 
(3) and (6) for the diffusion coefficient D are measures of the rate of energy diffusion, 
and Dq^, Vq^ are measures of the time scale for particle acceleration. Substituting the 
result (11) for kmin into equations (5) and (8), we find that Dq and Vq are given by 

for 2 < g < 4; 

V, = 'J^ (^) .h,aMi = ^^ f ^) '" ■hvil,Ra^'\ (14) 

8 \meJ 8 \"^p/ 

for 1 < g < 2. 

Corresponding to the Kolmogorov turbulent spectrum (g = 5/3), the function Jy/ is of 
order unity. As an idealized assumption, we regard the Kolmogorov spectrum as the 
representative spectrum over the range 1 < g < 2, and we henceforth set g = 5/3 and 
Jvr = 1 in (14). 

Since, from (10), the parameter a is inversely proportional to the particle number 
density A'^o, it follows from (13) and (14) that Dq and "Dq increase as Nq decreases. This 
agrees with the conclusions of Summers et al. [1998] who found by constructing resonant 
diffusion curves in velocity space that energy diffusion becomes more pronounced with 
increasing a (or decreasing A'"o). As expected, the values of Dq and Vq also increase as 
the turbulent spectral energy density ratio R increases. Specifically, we find from (13) 
and (14) that -Do and Vq depend on the plasma parameters A'o and i?o, and the wave 
amplitude Ai? as follows: 

D^o^B'i-''\/^Bf/Nt'^'^ (15) 

for 2 < g < 4; 

V,^Bl{^Bf/Nl'^ (16) 
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for 1 < g < 2. 

In this paper, we set A^o = 10 cm"'^ as the particle number density representative 
of the inner magnetosphere (3 < L < 9) outside the plasmasphere. It could be argued 
that such a value may be too high for the background plasma outside the plasmasphere. 
However, since from (15) and (16) it is clear that the acceleration process becomes 
more efficient as A^o decreases, we find it useful to adopt A'"o = 10 cm~^ as a generic 
conservative value. We comment further on this assumption below. We use the 
equatorial (dipole) magnetic field value Bq = 3.12 x 10^^/L^ T. Corresponding values 
of Bq and the above-defined parameters a and j3a at the locations L = 3, 4, ■ ■ -, 9 are 
given in Table 1. In Figure 0, we plot the energy diffusion parameter Dq (sec~^) as a Figure 
function of the spectral index q in the range 2 < g < 4 at each of the locations L = 3, 4, 
5. At each L-value, we calculate Dq for the specified wave amplitudes AB = 75 pT, 100 
pT, 300 pT, and 1 nT ( which correspond to the indicated values of R in the diagrams). 
Lines indicating the time scales Dq^ for particle acceleration corresponding to 1 hour, 
1/2 day, and 1 day are shown in each diagram. As can be observed from the curves in 
Figure 2, the value of Dq is particularly sensitive to the value of g as g approaches 2. 
In fact, formally from (13) we have the result Dq -^ oo, as g — *> 2, which is obviously 
undesirable physically, but which is a consequence of the quasi-linear diffusion formalism 
we have adopted in this paper. It is evident from Figure 2 that, at any given L-value, as 
the value of g decreases, the value of Dq increases, and hence the time scale for particle 
acceleration decreases. In addition, it can be observed that for a given value of g, as 
L decreases, the value of Dq likewise increases; this property also follows from relation 
(15). Thus, for g in the range 2 < g < 4, shorter acceleration times are favoured by 
smaller value of g, small values of L, and (of course) larger values of the wave amplitude 
AB. Corresponding to the Kolmogorov spectrum (g = 5/3), we show in Figure |^ the Figure ^ 
diffusion parameter T>q (sec~^) as a function of AB (pT) at the locations L = 3, 4, 5. 
Again, it is clear from the figure that shorter acceleration times are favoured by smaller 
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values of L and larger values of Ai?; this property similarly follows from (16). 



In Figure ^, we plot the diffusion coefficient D (sec~^) as a function of the particle Figure ^ 
kinetic energy E (MeV), as given by (3), (6), (13) and (14), for the fixed wave amplitude 
/\B = 1 nT, and A^o = 10 cm~^. The curves are constructed for values of the spectral 
index q = 5/3, 2.5, 3, 3.5, and 4, at each of the locations, L = 3, 4, 5. The diffusion 
coefficient D is clearly an increasing function of energy E. In general, for a given value 
of E, though not at all values, D can also be seen to increase as q decreases. 

According to the standard quasi-linear theory of resonant interaction of electrons 
with whistler-mode turbulence [e.g., Melrose, 1986], in order for electrons to resonate 
with whistlers, the condition, 

7/5 > {mp/m,y/^(3A (17) 

must be satisfied; P = v/c where v is the particle speed, c is the speed of light, 7 is the 
Lorentz factor, and jSa = va/c is the Alfven speed parameter defined above. Making use 
of relativistic relations given in (A6), we find that (17) can be expressed in the form 



E{E + 2) > {mp/me)Pl 



:i8) 



which, in turn, by using the parameter a defined in (10), can be reduced to 



E> Er 



(19) 



where Ec is the critical energy given by 



E, = (l + aY/'-l. 



(20) 



The value of the parameter a depends on the values of the particle number density 
A^o and magnetic field Bq. Values of the critical energy E^ are given in Table 1 at the 
locations L = 3, 4, • ■ -, 9, for Nq = 10 cm~^; the values for L = 3, 4, 5 correspond 
respectively to the energy cut-off values in the upper, middle, and lower diagrams in 
Figure 4. 



Table 1 
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3. NUMERICAL RESULTS 

Prior to consideration of the solution of the kinetic equation (1) for the electron 
energy distribution function f{E,t), we must specify the source function S{E,t) which 
represents storm-produced seed electrons. We shall assume that the source function can 
be represented by a standard relativistic Maxwellian distribution, namely, 

S = So[fi/KM]{E + 1)[E{E + 2)]V2e-Mi^+i) (21) 

where 

/i = rrieCyikBT,) (22) 

and Ts represents the temperature of the distribution; K2{fi) is a modified Bessel 
function of the second kind of argument /i, and ks is Boltzmann's constant. It can be 
shown that / SdE = Sq, so that the parameter 5*0 represents the total number of source 
electrons per unit volume per unit time. 

It is clearly important both to calibrate and test our model, as far as is possible 
at the present time, by making use of available observational data. The study by 
Cayton et al. [1989] appears best-suited to these purposes. Cayton et al. [1989] derived 
energy distribution functions from energetic (30 - 2000 keV) electron fluxes observed 
simultaneously by three satellites in geosynchronous orbit throughout the year 1986. 
It was found that the energetic electron population can be resolved into two distinct 
relativistic Maxwellian components, each parameterized by a temperature and a density: 
a lower energy (30 - 300 keV) "soft" electron distribution with temperature Tg ^ 25 
keV and number density A^g ^ 5 x 10^^ cm^^; and a higher energy (300 - 2000 keV) 
"hard" electron distribution with temperature Th ~ 200 keV and number density 
Nh ^ 10~^cm^'^. The "soft" component is charaterized by intense substorm-related 
injections and by strong temporal variations. Accordingly, and in agreement with a 
suggested interpretation by Cayton et al. [1989], we shall regard this "soft" component 
as comprising the electron seed population. Thus, we shall identify the temperature Tg 
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associated with the Maxwelhan source distribution (21) - (22) as the aforementioned 
temperature of the "soft" electron component. Cayton et al. [1989] found that the value 
of Th shows little change on the substorm (hourly) time scale, while Nh decreases during 
substorms. We shall regard the "hard" electron distribution as precisely the highly 
energetic ("killer") electron distribution that we are trying to model as a (steady-state) 
solution of the kinetic equation (1) with the steady Maxwellian source (21) - (22). 

We solve equation (1) for the energetic electron distribution f{E,t) by the Crank- 
Nicholson implicit differencing scheme. The method is well suited to time- dependent 
Fokker-Planck equations, and we refer the reader to Hamilton et al. [1990] and Park 
and Petrosian [1996] for full details. Since we are concerned with the generation of 
a highly-energetic electron distribution, we assume that there are no such energetic 
particles initially, i.e., 

/(E,0) = 0, E>Es, (23) 

where Es = l//i is the thermal energy associated with the source distribution (21) - 
(22). We further assume that, subject to continuous injection of the seed electrons 
(given by (21) - (22)), the evolving distribution maintains a maximum at E = Es for all 
time, i.e., we take the inner boundary condition as 

m^ = 0. E = E,. (24) 

dE ' " ^ ' 

Finally, for the outer boundary condition we require that the distribution function tend 
to zero for large values of E for all time, so we set 

f{E,t) = 0, E>Eo (25) 

where Eq is a specified upper value of E (in practice, we fix i?o = 2 x 10^ MeV). 
Having constructed the evolving electron distribution subject to the above conditions, we 
thereby obtain the resulting steady-state distribution f{E) which we fit to a relativistic 



Maxwellian distribution, i.e., we carry out the linear fit. 



16 



logio \f{E)/{{E + l)[E{E + 2)] V2}] ^a + hE 



(26) 



with 



and 



a = loE 
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iV^Ae-Vi^2(A) 



-Alogioe, 



(27) 



A = mecyiksn), (28) 

where Th and A^/i represent respectively the temperature and number density of the 
steady-state distribution (to be compared with the above values associated with the 
"hard" electron distribution); and K2{X) is a modified Bessel function of the second kind 
of argument A. For a given steady-state solution, the parameters a and b are determined 
by a linear regression comprising a minimization of a chi-square goodness-of-fit merit 
function [Press et al, 1992]. Having thus obtained values for a and b, we then calculate 
Th and A^^ from (27) and (28). 

Representative numerical solutions of the model presented in this paper are shown 
in Figures ^, |^, and |^ with corresponding results given respectively in Tables 2, 3, 
and 4. In all cases we set the background number density A^o = 10 cm~^, and the 
source electron temperature T^ = 25 keV (giving /i = 20.44 from (22)), the latter value 
being equal to the estimate by Cayton et al. [1989] for the temperature of the "soft" 
electron distribution. The scheme and rationale for setting the remaining parameters 
is as follows. Firstly, we set L, which fixes the value of the background magnetic 
field Bq. Secondly, we set the average wave- amplitude AB and the spectral index q; 
these values are chosen to yield a value of the diffusion parameter Dq (or Vq) that is 
expected to produce a steady-state (equilibrium) distribution function after several 
days of source injection. Whistler-mode "chorus" wave amplitudes have been reported 
in the range 1 - 100 pT [Burtis and Helliwell, 1975], with Parrot and Gaye [1994] 
finding that during more intense periods of magnetic activity wave amplitudes can 
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approach AB = 1 iiT. Amplitudes of whistlers associated mainly with hiss, with values 
of 100 pT or more, have also been reported by Smith et al. [1974] during a typical 
storm recovery phase. In Table 1, corresponding to a background number density of 
A^o = 10 cm~^, we present the wave amplitudes AB (pT) expected to yield a high 
energy "hard" electron distribution after a few days of seed electron injection as a result 
of substorm activity. These required amplitudes depend on the value of L and q (as well 
as Nq). For q = 5/3, and for the inner region 3 < L < 5, the values of AB are in the 
range 75 - 400 pT, which are realistic though in the higher range of observations. As 
pointed out in Section 2, the process of gyroresonant stochastic acceleration becomes 
more efficient with decreasing background plasma density. While Nq = 10 cm^'^ can 
be regarded as a representative value for the background plasma density outside the 
plasmasphere in certain conditions, it is also true that at other times Nq = 1 cm~^ is a 
more representative value. Using (15), we calculate that if we set A"o = 1 cm~'^, then, 
for 3 < L < 5, the required Ai?-values are in the range 13 - 72 pT, if g = 5/3, and in 
the range 35 - 316 pT, if 2.5 < q < 3. Having set values for No,Ts,L,q, and AB, we 
next specify the mean particle escape time Tesc- In fact, since a value for Tesc is not 
precisely known, we treat Tesc as an adjustable parameter and run the cases l/(-Do^esc) 
(or l/(X'oTesc)) = 0,1,2,5,10. In Figure 5 (left) we show steady-state solutions for 
the electron energy distribution function f{E) for the case L = 3, q = 5/3, AB = 75 
pT. The source strength Sq has been chosen so as to produce a model solution that 
best agrees with the "hard" electron distribution of Cayton et al. [1989]. In order to 
achieve this, for each steady-state solution the linear fit (26) - (28) to a relativistic 
Maxwellian distribution is carried out. The corresponding results are shown in Figure 
^ (right) and Table 2. It is found that the temperature Th associated with a particular 
steady-state solution, as given by the parameter b (or the slope of the constructed line), 
is largely determined by the value of Tesc, while the number density Nh, which is then 
given by the parameter a (or the vertical intercept of the line) is largely determined by 



the value of Sq. The results in Table 2, for which 5*0 = 1.5 x 10~^ cm~^sec~^, indicate 
a best agreement with the empirically-derived values of Nh ~ 10~^ cm~'^ and Th = 200 
keV, when Tesc ~ 1/2 day, the corresponding time for the formation of the steady-state 
solution being Teq ~ 4 days. A precise, physically representative value for the mean 
particle escape time T^sc is difficult to determine a priori since T^sc relates to scattering 
losses of electrons due to both whistler-mode and EMIC waves. However, on the basis 
of estimates of timescales for strong diffusion scattering loss, it appears that Tgsc is of 
the order of hours and so T^sc ~ 1/2 day is not an unreasonable value. We relate Teq to 
the time taken after the initiation of a storm for fluxes of relativistic electrons to peak 
(see Section 1), which is observed to be several days. Thus, we favour solutions of the 
present model for which Teq =1-5 days, with Teq ~ 4 days possibly the optimal 
value. 

In Figure ^ and Table 3 we show the corresponding results for the case A^o = 10 
cm~^, L = 6.6, q = 2.5, AB = 800 pT, while in Figure ^ and Table 4 we show the results 
for the case A'"o = 10 cm~^, L = 5, q = 3, AB = 1 nT. As can be seen from the tables, 
for both these cases best agreement between the solutions and the "hard" electron 
distribution of Cayton et al. [1989] occurs when Tesc ~ 1/2 day, and corresponds to a 
formation time Teq = 3-5 days. If the background number density is taken to be 
Aq = 1 cm~^, the cases shown in Figures 6 and 7 correspond respectively to values for 
the wave amplitude AB of 190 pT and 316 pT. Figures 6 and 7 correspond to cases of 
intense substorm activity during the storm recovery phase. 

Taking into account the value of the wave amplitudes given in Table 2 corresponding 
to A'o = 10 cm~^, and their converted values for the case Aq = 1 cm~^, we re- iterate 
that the model solutions imply in particular that for a Kolmogorov turbulent wave 
spectrum, sustained whistler amplitudes in the physically reahstic range 13 - 72 pT can 
generate a typical high energy "hard" electron distribution in the inner magnetosphere 
3 < L < 5 within one or two days. It should also be noted that the model calculations 
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show that the acceleration mechanism considered in this paper is not effective in the 
region 7 < L < 9 since the necessary values of the whistler amplitude would be too 
high (typically in excess of 1 nT). Thus, the model formulated herein has been shown 
to be a viable mechanism for accelerating electrons exactly in the inner region of the 
magnetosphere where the peak in electron phase space density of the highly energetic 
electrons is observed to occur [e.g., Selesnick and Blake, 1997b]. 

A requirement of the model presented here is enhanced whistler-mode chorus 
lasting for a period of at least one or two days. Geomagnetic conditions during which 
such a requirement is particularly well satisfied occur during the descending phase of the 
solar cycle when the Earth's magnetosphere can be impacted by a high-speed solar wind 
stream following a magnetic field build-up known as a Corotating Interaction Region 
(CIR). CIRs cause small and moderate magnetic storms, but not major storms. Since 
the Earth can be embedded in the associated high-speed stream for days to weeks there 
are substorms for days to weeks [Tsurutani et al, 1995; Kamide et al, 1998]. Thus, 
during this long-lasting recovery phase of the magnetic storm, there will be continuously 
enhanced wave activity, in the form of both whistler-mode chorus and EMIC waves, 
to drive the acceleration mechanism presented herein to generate the high energy (>1 
MeV) electrons. 

4. CONCLUSIONS 

In this paper, by means of quasi-linear theory and a test-particle approach, we have 
formulated the model kinetic equation (1) in which the acceleration mechanism is due to 
gyroresonant interaction between electrons and whistler-mode turbulence, corresponding 
to parallel wave propagation. The essential purpose of the study has been to apply 
equation (1) to the Earth's inner magnetosphere in order to test the hypothesis that 
storm-enhanced whistler-mode chorus can accelerate lower-energy substorm-produced 
seed electrons to relativistic (> 1 MeV) energies over a period of a few days. Our 



20 



conclusions are as follows: 

1. Based on the model calculations in this paper, it is entirely possible for enhanced 
whistler-mode chorus to generate the observed increases in relativistic (> 1 
MeV) "killer" electrons during the storm recovery after a period of several days, 
so long as the waves are sufficiently strong. If Nq = 10 cm^'^ is taken to be 
the background plasma number density outside the plasmasphere, the typical 
average wave-amplitudes required for a Kolmogorov spectrum are in the range 
AB = 75 - 400 pT, dependent on the location L. If A'^o = 1 cm~^, the required 
wave-amplitudes are in range 13 - 72 pT. 

2. Energetic electron spectra of the model solutions do not follow a simple power law 
in energy. For certain sets of parameters, we find that solutions can be well fitted 
to the relativistic Maxwellian distribution empirically constructed by Cayton et 
al. [1989] to represent the higher energy (300 keV - 2 MeV) "hard" electron 
population at geosynchronous orbit. We note the recent analysis by Freeman et al. 
[1998] of the November 3-4, 1993 storm, in which electrons from about 100 keV to 
1.5 MeV were characterized by a power law spectrum. Evidently, optimal fitting of 
empirical electron spectra to power-law, Maxwellian, or other types of distribution 
can depend on the energy range prescribed and the event under consideration. In 
connection with electron power-law energy spectra. Ma and Summers [1998] have 
shown that such spectra can be produced by whistler-mode turbulence, although 
it is questionable whether the necessary conditions established in their theoretical 
study can be satisfied in the Earth's magnetosphere. 

3. It is unlikely that any single physical mechanism of electron acceleration can 
fully account for relativistic electron enhancements occurring during the recovery 
phase of magnetic storms, not least because various types of energetic electron 
event have been observed [e.g.. Baker et a/., 1997, 1998; Reeves, 1998b; Reeves et 
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ai, 1998]. Rapid energetic electron flux enhancements taking place over minutes 
have been associated with inductive electric fields [e.g., Li et ai, 1993], while 
enhancements occurring over tens of minutes or a few hours have been linked 
to ULF pulsations [e.g., Rostoker et ai, 1998; Liu et al, 1999]. The gradual 
acceleration process (occurring over a few days) formulated in this paper is not 
intended to apply to such energetic electron events which typically result from 
major storms. However, small and moderate magnetic storms associated with 
Corotating Interaction Regions (CIRs) characteristically have long recovery phases 
and attendant substorms for days to weeks [Tsurutani et al., 1995; Kamide et al., 
1998]. Since these substorms produce enhanced whistler-mode chorus (and EMIC 
waves) over (at least) several days, the necessary conditions for the effectiveness 
of the mechanism presented in this paper are satisfied. Hence, for these types of 
storm, and possibly others, when average wave-amplitudes are sufficiently large, 
the present study shows that, in conjunction with pitch-angle scattering by EMIC 
waves, the mechanism of stochastic acceleration by whistler-mode turbulence is a 
serious candidate for explaining the generation of "killer electrons" . 

Appendix A: Derivation of the Fokker-Planck equation (1) 

We consider energetic charged particles in a uniform magnetic field, with 
superimposed small-amplitude plasma waves of a given mode. The equation for the 
evolution of the particle distribution function 0(p, t, jj) due to gyroresonant interactions 
with the waves is 

dt p^ dp \ ^^ dp J p2 Qp y ^^ diij dfi \ ^^ dp J 

Equation ([Al|), called a kinetic or diffusion or Fokker-Planck equation, is derived by 
expanding a collisionless Boltzmann equation for the particle distribution function 
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to second order in perturbed quantities, and ensemble-averaging over the statistical 
properties of the plasma waves in accordance with quasi-linear theory. Among the 
early authors to carry out this procedure were Kennel and Engelmann [1966], Hall and 
Sturrock [1967], and Lerche [1968]; see also Melrose [1980], Schlickeiser [1989], and 
Steinacker and Miller [1992], and references therein. In ([Al]), p is the relativistic unit 
momentum given by p = 'jv/c, where v is the particle speed, and 7 = (1 — v'^/c^)~^^'^ is 
the Lorentz factor, with c the speed of light; t is time; fi is the cosine of the pitch angle; 
Sl is an energy loss term due to processes not directly associated with gyro-resonant 
wave-particle interactions; and Q{p,t) is a source term. The Fokker-Planck or diffusion 
coefficients Dpp, Dp^, D^p, and D^^ depend on the properties of the wave turbulence, 
viz., the wave mode and polarization, the angle of wave propagation to the ambient 
magnetic field, and the power spectrum, including the ratio of the turbulent wave 
energy to the background magnetic energy. These coefficients have been given both 
in general form, and specific form, for various particular wave modes, by a number 
of authors, e.g., Melrose [1980], Schlickeiser [1989], Steinacker and Miller [1992], and 
Hamilton and Petrosian [1992]. It is not necessary here to derive equation (|A1|) , which 
requires considerable algebra, or to provide general expressions for the coefficients 
Dpp, Dp^ = D^p, and D^^. We shall assume that the rate of pitch-angle scattering is 
much larger than the rate of energy diffusion (and the rate of the particle escape from 
the system). Such an assumption is reasonable based on an analysis of time scales 
associated with resonant interaction of electrons with whistler-mode waves, e.g., see 
Melrose [1980] and the discussion by Hamilton and Petrosian [1992]. Equivalently, 
defining the time scales T^^ = D~^, Tfj,p = pD^^, Tpp = p'^D~p, and the escape time Tesc, 
we assume that T^^ <^ Tpp, T^^ <^ T^p, and T^^ <^ Tesc- Then the particle distribution 
function can be assumed to be isotropic, and the pitch angle can be eliminated from the 
problem by integrating ( [Al] ) with respect to /x [e.g., see Schlickeiser, 1989; Steinacker 



and Miller, 1992]. Writing 
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F{p,t)= I ^<j){p,t,fi)dfi, (A2) 



and representing the scattering loss of particles by pitch angle diffusion by means of a 
loss term —F{p,t)/Tesc, the equation (|7n|) thus becomes 



9F{p,t) ^ I d (^.j^ ,. dF{p,t) 
dt p"^ dp \ ^ dp 



)+^|(^'^^^^)^(^'^ 



^^^'^^+iQ(p,t), (A3) 



T 2 

-'■ esc ^ 

where the momentum diffusion coefficient Dp{p) has been formed by averaging with 
respect to /x. 

We now change the momentum variable p to the kinetic energy variable E = j — 1 
in equation ([A3|) . We write 

f{E, t)dE = 4Vf(p, t)dp, (A4) 

^ = ^ A, (A5) 

dp dp dE ' 

and make note of the following simple relativistic relations: 

p = -fP,P = v/c,p^ = E{E + 2),j = {1 + p^y/\ 
pdp ={E + l)dE, (3dp = dE,(3= [E{E + 2)Y/\E + 1)-^ (A6) 



Then, after straightforward manipulation, equation ([A3|) can be expressed in the form 



I {fiE, t)) = A_ [DiE)f{E, t)] - A [(^AiE)- I E, |) f{E, t)] - l^ + 5(E, t). 



- esc 



(A7) 
where 

D{E) = P'Dpip), 

\El\ = PSl{p), 
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S{E,t) = ^Q{p,t). (A8) 

The form of equation ( |A7|) is actually the Fokker-Planck form of equation for a particle 
distribution function as originally presented by Chandrasekhar [1943] for particles in 
stochastic motion. Stochastic acceleration of electrons in solar flares has been treated 
using different versions of ( |A7| ), e.g., see Ramaty [1979], Petrosian [1994], and Park et 



al. [1997]. 
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Figure 1. (Left) Schematic view in the magnetic equatorial plane of the approximately 
circular (projected) drift path of relativistic electrons in the inner magnetosphere. During 
storms these energetic electrons drift (eastward) through regions of enhanced whistler- 
mode chorus and enhanced electromagnetic ion cyclotron (EMIC) waves. (Right) Repre- 
sentation of the gyration about magnetic field lines and the bounce motion of energetic 
electrons as they execute the approximately circular drift path shown in (a). 
Figure 2. Diffusion coefficient Dq given by (13) as a function of the turbulence spectral 
index q, for 2 < g < 4. The upper, middle, and lower diagrams correspond respectively to 
the locations L = 3,4, 5. The background particle number density A^o = 10 cm~^. In each 
diagram, curves are shown corresponding to the four indicated values of the wave power 
R (given by (9)) which correspond to the respective average wave amplitudes AB =75 
pT, 100 pT, 300 pT, 1 nT. 

Figure 3. Diffusion coefficient Vq given by (14) as a function of the average wave 
amplitude AB (pT), at each of the locations L = 3,4,5. The turbulence spectral index 
q = 5/3, the parameter J^ = 1, and the background particle number density Nq = 10 

cm~^. 

Figure 4. Diffusion coefficient D given by (3), (6), (13), (14) as a function of the particle 

kinetic energy E, for the average wave amplitude AB = 1 nT, and the indicated values 

of the spectral index q. The upper, middle, and lower diagrams correspond respectively 

to the locations L = 3, 4, 5. The background particle number density A^o = 10 cm~^. 
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Figure 5. (Left) Steady state solutions f{E) to tlie kinetic equation (1) for the electron 
energy distribution function, for the indicated values of l/iVoTesc), corresponding to 
different mean particle escape times. The diffusion coefficient and systematic acceleration 
rate are given by (6) and (7), with the diffusion parameter Vq defined by (f4); q = 5/3, 
Jw = 1, L = 3, AB = 75 pT, A^'o = fO cm~^ and Vq = 7.8 x fO^^ sec^^ The particle 
source function is given by (21), with 5*0 = 1.5 x 10~^ cm~^ sec~^ and /x = 20.44. (Right) 
Corresponding re-scaled plots of the solution curves on the left, for comparison with 
relativistic Maxwellian energy distribution functions. The dashed lines represent best 
fits with Maxwellian distributions in accordance with the results given in Table 2. 
Figure 6. (Left) Steady state solutions f{E) to the kinetic equation (1) for the electron 
energy distribution function, for the indicated values of l/(Do^esc)) corresponding to 
different mean particle escape times. The diffusion coefficient and systematic acceleration 
rate are given by (3) and (4), with the diffusion parameter Dq defined by (13); q = 2.5, 
L = 6.6, AB = 800 pT, A^o = 10 cm^^ and Do = 7.2 x 10"^ sec^^ The particle source 
function is given by (21), with 5*0 = 1.4 x 10^^ cm^'^ sec"^ and fi = 20.44. (Right) 
Corresponding re-scaled plots of the solution curves on the left, for comparison with 
relativistic Maxwellian energy distribution functions. The dashed lines represent best 
fits with Maxwellian distributions in accordance with the results given in Table 3. 
Figure 7. (Left) As in Figure 6 (Left), but for the parameters q = 3, L = 5, AB = 1 
nT, A'o = 10 cm^^, and Dq = 6.5 x 10^^ sec^^. The particle source function is given by 
(21), with 5*0 = 1.9 X 10^^ cm^^ sec^^ and fi = 20.44. (Right) Corresponding re-scaled 
plots of the solution curves on the left, for comparison with relativistic Maxwellian energy 
distribution functions. The dashed lines represent best fits with Maxwellian distributions 
in accordance with the results given in Table 4. 
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Table 1. Values of the magnetic field Bq (10^^ T), the parameters 
a and j3a given by (10), and the critical energy E^ (keV) given by 
(20), corresponding to the locations L = 3, 4, ■ ■ -, 9; the background 
particle number density is A^o = 10 cm^^. Also given, corresponding 
to the indicated values of the spectral index q, are typical values of 
the average wave amplitude AB (pT) and the associated values of the 
diffusion parameters "Dq and Dq in sec~^(xlO~^) required to produce a 
high energy "hard" electron distribution after several days of sub-storm 
particle injection. 





Bo 


a 


(3a 


Ec 


g = 5/3 


g = 2.5 


q = 3.0 


L 


AB Do 


AB Do 


AB Do 



3 11.6 1.31 0.027 267 75 7.8 150 8.6 500 7.6 

4 4.85 0.23 0.011 57 200 9.8 300 9.3 700 6.2 

5 2.50 0.061 0.006 15 400 10 400 6.1 1000 6.5 

6 1.44 0.019 0.0032 4.8 600 7.3 600 5.7 1200 5.3 

7 0.91 0.0081 0.0021 2.1 900 7.0 800 5.4 1500 5.4 

8 0.61 0.0036 0.0014 0.9 1400 7.5 1200 6.6 1800 5.1 

9 0.43 0.0018 0.001 0.5 2000 7.6 1500 6.1 2200 5.4 
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Table 2. Results associated with the steady-state solutions shown in Figure 5. Each line of 
the table corresponds to the particular value of T^sc (the mean particle escape time) indicated. 
Each solution is fitted to a relativistic Maxwellian distribution by means of the linear fit (26) in 
which the parameters a and b yield values for the number density A'^^ and temperature Th of the 
distribution; x^ measures the goodness-of-fit. The time taken for the steady-state (equilibrium) 
distribution to form is Teq- 
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Table 3. As for Table 2, except the results are associated with the steady-state solutions in 
Figure 6. 
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Table 4. As for Table 2, except the results are associated with the steady-state solutions in 
Figure 7. 
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